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HIGH ORDER NUMERICAL SCHEMES FOR SECOND-ORDER 
FBSDES WITH APPLICATIONS TO STOCHASTIC OPTIMAL 

CONTROL * 

TAO KONGt, WEIDONG ZHAO*, AND TAO ZHOU § 

Abstract. This is one of our series papers on multistep schemes for solving forward backward 
stochastic differential equations (FBSDEs) and related problems. Here we extend (with non-trivial 
updates) our multistep schemes in [W. Zhao, Y. Fu and T. Zhou, SIAM J. Sci. Comput., 36 (2014), 
pp. A1731-A1751.] to solve the second order FBSDEs (2FBSDEs). The key feature of the multistep 
schemes is that the Euler method is used to discrete the forward SDE, which dramatically reduces 
the entire computational complexity. Moreover, it is shown that the usual quantities of interest (e.g., 
the solution tuple (Yt , Zt, At , Tt) in the 2FBSDEs) are still of high order accuracy. Several numerical 
examples are given to show the effective of the proposed numerical schemes. Applications of our 
numerical schemes for stochastic optimal control problems are also presented. 

Key words. Multistep schemes, the Euler method, second-order forward backward stochastic 
differential equations, stochastic optimal control 
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1. Introduction. This work is concerned with numerical methods for the fol¬ 
lowing coupled second order forward backward stochastic differential equations (2FB- 
SDEs) which is defined on the filtered probability space (12, J^", F, P) : 

X t = x+ f 6(s,0 s )ds+ f a(s, 0 s )dPF s , 

Jo Jo 

f{s,e 3 )ds-J t z s dw„ s £ [0,T], 

Z s = Zq + f A s ds + f r g dW s , 

Jo Jo 

where 0 t = Z t , A t ,T t ) £ K m x R x x S d is the unknown, (f2 ,&,P) is 

the given probability space, T > 0 is the deterministic terminal time, [Wt}t^\o,T] is 
a d-dimensional Brownian Motion defined on (12,^", P) with the natural filtration 
F = {J^t}o<t<T and all P-null sets in J^o, x £ is the initial condition of the 
forward SDE, S d is the set of all d x d real-valued symmetric matrices, and 

b : [0, T] x R m x I x R d x S d M m , a : [0, r]xK m xRxK < ‘x5 i -> R mxd 

are referred to the drift and diffusion coefficients of the forward SDE, respextively. 
While 


( 1 . 1 ) 


Y t = g(X T ) + 


f : [0, T] x x I x R d x S d —!> M, g : R m -A K 
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are referred to the generator and the terminal condition of the backward SDE, respec- 
tively. The two stochastic integrals with respect to the Brownian Motion {W t }t£[o.T] 
are of the Ito type. A 5-tuple (X t , Y t , Z t , T t , A t ) is called an A 2 -adapted solution of 
the 2FBSDEs (1.1) if it is ^-adapted, square integrable, and satisfies (1.1). More¬ 
over, the 2FBSDE (1.1) is called decoupled if b and er are independent of Y t , Z t A t 
and T t . 

The 2FBSDEs (1.1), as an extension of the BSDEs in the linear [2] or nonlinear 
cases [19], was first introduced by the cornerstone work of Cheridito, Soner, Touzi 
and Victoir [5], with a slightly different (yet equivalent) formula, and it was further 
investigated by Soner, Touzi and Zhang in [22]. The main motivation there is to give 
a precise connection between the 2FBSDEs and fully non-linear PDEs, in particular 
the Hamilton-Jacobi-Bcllman equations and the Bcllman-Isaacs equations which are 
widely used in stochastic control and in stochastic differential games. Such a con¬ 
nection leads to interesting stochastic representation results for fully nonlinear PDEs, 
generalizing the original (nonlinear) Feynman-Kac representations of linear and semi- 
linear parabolic PDEs, see e.g. [11,16,20,21] and references therein. 

From the numerical point of view, one can adopt these connections between PDEs 
and FBSDEs to design the so called probabilistic numerical methods for PDEs, by 
solving the equivalent FBSDEs (or 2FBSDEs). While there are a lot of work dealing 
with numerical schemes for BSDEs [1,3,4,6,13,24,25,27], however, there are only a 
few work on numerical methods for FBSDEs [7,9,14,15,17,23,26,28] and fully non¬ 
linear PDEs [8,10]. Some of the above works are designed with high order accuracy 
that can however only be used to deal with low dimensional FBSDEs. While some 
of them are low order numerical methods that are suitable for solving high dimen¬ 
sional problems. In particular, we mention the work [10], where a numerical example 
for a 12-dimensional coupled FBSDE is reported, and it is shown by numerical test 
that the numerical method converges with order 1. Also, in [9], multistep schemes 
were proposed to solve multi-dimensional FBSDEs by using the spares grid interpo¬ 
lation technique, and several multi-dimensional examples with dimension up to 6 are 
presented, and high order convergence rates up to 3 were obtained. 

To the best of our knowledge, there is no related studies for high order numerical 
methods for 2FBSDEs. The main purpose in this work is to extend our multistep 
schemes in [26] (which is original designed for solving FBSDEs) to the use of solving 
2FBSDEs. The key feature of the multistep schemes is that the Euler method is used 
to discrete the forward SDE, which dramatically reduces the entire computational 
complexity. Moreover, it is shown that the quantities of interest (e.g., the solution tu¬ 
ple (Y t , Z t , r t , A t ) in the 2FBSDEs) are still of high order accuracy. Several numerical 
examples are given to show the effective of the proposed numerical schemes. With 
particular interests, applications of the proposed numerical schemes to stochastic op¬ 
timal control problems are also presented and corresponding numerical examples are 
also given. It is worth to note that the multistep schemes proposed in this paper can 
also be used to solve a large class of fully nonlinear PDEs and we left this for our 
future studies. 

The rest of the paper is organized as follows. In Section 2, we present some 
preliminaries. Then in Section 3, We shall discuss the multistep schemes for solving 
decoupled 2FBSDEs, and this is followed by extensions to coupled 2FBSDEs in Section 
4. In Section 5, we shall present several numerical examples, and applications of our 
numerical methods for stochastic optimal control problems will also be presented. We 
finally give some concluding remarks in Section 6. 
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2. Preliminaries. We first introduce some notations that will be used in this 
paper. For x £ R d , B £ S d , we set |x| := y/x'f-\ - \-x 2 d , \B\ := Bfj- We 

denote by C k the set of functions g(x) : K —x R with uniformly bounded derivatives 
up to order k, and by C kl,k2 the set of functions f(t,x) : [0, T] xl" 4 R m with 
continuous partial derivatives up to order k\ w.r.t. t £ R and up to order &2 w.r.t. 
x £ R m . 

2.1. Diffusion process and its generator. A stochastic process X t is called 
a diffusion process starting at xo and at the time to if it satisfies the following SDE 

(2.1) X t =x 0 + [ b(s,X s )ds+ [ a{s,X s )dW s , t £ [to,T], 

Jto Jta 

where b s = b(s,X s ) and a s = a(s,X s ) are measurable functions satisfying 

|6(s,x)| + |ct(s,x)| < C(1 + |x|), x £ R m ,s £ [t 0 ,T], 

I b{t,x) - b(t,y) | + |cr(£, x) - <r{t,y)\ <L\x-y\, x,y £ K m ,s £ [t 0 ,T], 


It is well known that under conditions (2.2), the SDE (2.1) admits a unique solution. 
It is worth to note that, by the Markov property of the diffusion process, we have 
Ef [V s ] = E [X s | X t = x\ ,Vt < s. 

Let Xt be the solution of (2.1). Then, for any given measurable function g : 
[0, T] x R m —>• R, g(t, X t ) becomes a stochastic process. Moreover, we give the follow¬ 
ing definition 

Definition 2.1. Let X s be a diffusion process in that satisfies (2.1). The 
generator A of X s on g is defined by 


(2.3) 


Ag(t, x) = lim 

S^.t 


[flC^ATs)] — git, x) 
s — t 


x £ R”. 


Concerning the generator A, the following result holds [18]: 

Theorem 1. Let X s be the diffusion process defined by the SDE (2.1). If f £ 
C 1,2 ([0,T] x R m ), then we have 


(2.4) Af(t, x) = Cf{t, x), Af(t, X t ) = Cf(t, X t ), 


where the operator C is defined by 

(2.5) £cj>(t,x) := <j>t{t,x) + S7 x cj)(t, x)b(t, x) + ^tr (a(t, x)a T (t, x)X 2 x cj)(t, x)) . 

with X x (j) = ( d Xl (f i,..., d Xm (f>), and X x (f> being the Hessian matrix of <j> with respect to 
the spatial variable x. 

Note that Af(t,X t ) £ T t is a stochastic process. Furthermore, by using together 
the Ito’s formula, Theorem 1 and the tower rule of conditional expectations, we have 
the following theorem. 

Theorem 2. Let to < t be a fixed time, and xo £ R m be a fixed space point. If 
f £ C 1,2 ([0, T] x R m ) and E“° [\Af(t,X t )\] < +oo, we have 

(2.6) ^d =E %[Af(t,X t )], t>to. 
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Moreover, the following identity holds 


(2.7) 


dE*° [f(t,X t )} 


dE^ o ° [f(t,X t )] 


d t 


d t 


where Xt is a diffusion process satisfying 


( 2 . 8 ) 


X t — x o + 





with b : [0, T\ x R m —> R. m , a : [0, T] x R m —> R mxd being smooth functions satisfying 


b(t 0 , x 0 ) = b(t Q , x 0 ), cr{t 0 , x 0 ) = a(t Q , x 0 ). 


Note that by choosing different b and a, the identity (2.7) yields different ways for 
approximating dE^ o [f(t,X t )] /dt\ t=t . The computational complexity can be signifi¬ 
cantly reduced if appropriate choices are made. 

2.2. Solution regularity and representation theory of 2FBSDEs. We 

now consider the following decoupled 2FBSDEs 



(2.9) 


with a given terminal condition Yj- = g{X t). Under some standard assumptions (for 
details, please refer to [5,22]), the following results are shown [5]: 

Theorem 2.2. Let u = u(t,x ) be the solution of the following fully nonlinear 


PDE 


£u + f(t,x, u, V x ua, V x (X7 x ua)cr) = 0 


( 2 . 10 ) 


with the terminal condition u(T 7 x) = g{ x), and let (X t7 Y tl Z t ,T tl A t ) be the solution 
of the 2FBSDE (2.9). Then we have 


Y t = u(t,X t ), Z t = (V x ua)(t,X t ), 

T t = (V x (X x ua)a)(t,X t ), A t = (C(W x ua))(t,X t ). 


where the associate operator C is defined by (2.5). 

The above theorem provides a stochastic representation for solutions of fully non¬ 
linear parabolic PDEs, generalizing the pioneer work on Feynman-Kac representations 
of linear and semi-linear parabolic PDEs [11,16,20,21], 

2.3. Derivative approximation. We now recall some basic results for numer¬ 
ical approximation of derivatives, and these results will play an important role for 
designing our high order numerical methods for 2FBSDEs. 

Let u(t) £ Cjf +1 and U £ R (0 < * < k) satisfying to < ti < • ■ • < tfc, where k is 
a positive integer. Let Atoy = t» — to, * = 0,1,..., k. Then by Taylor’s expansion, for 
each tj, i = 0, 1 ,..., k, we have 



k 
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Then, we can deduce 


k 

k k _.j / k 

Y a k,Mti ) = Y —— j\ -+ ° (5Z a kA At o,i) k+1 

i —0 j —0 \i—0 

where otk,i, i = 0,1,..., k, are real nubmers. By choosing a^,i (i = 0,1,..., fc) as 



( 2 . 11 ) 


Yh a k,i( A to,iY 

i=0 _ 

j! 


1, 3 = 1 
o, 3*1 


One obtains a high order approximation 

k 


( 2 . 12 ) 


du 


^7 (to) = Y a k,i u (U ) +i?D, 


i=0 


where = 0(Y2i=o a fc,i(Ato,i) fc+1 )- In particular, when Af 0 ,i = iAt, we get from 
(2.11) the following linear system for c^yAt, 


(2-13) ^> fe ,At] = r’ 3 

ti l°> 3*1 

Note that the above system can be solved easily. We list akA At (i = 0,1 of 

the system (2.13) for fc = 1, 2,..., 6 in the following table. 


Table 2.1 


Oik 

,*A t 

i = 0 

i = 1 

i = 2 

i = 3 

2 = 4 

2 = 5 

2 = 6 

k 

= 1 

-1 

i 






k 

= 2 

3 

2 

2 

1 

2 





k 

= 3 

11 

6 

3 

3 

2 

1 

3 




k 

= 4 

25 

12 

4 

-3 

4 

3 

1 

4 



k 

= 5 

137 

60 

5 

-5 

10 

3 

5 

4 

1 


k 

= 6 

49 

6 

15 

20 

15 

6 

1 


20 

2 

3 

4 

5 

6 


The above approximation schemes are very popular in numerical ODE community, 
and it is known that such an approximation scheme is unstable for k > 7, and this is 
why we have only listed the values of afc^Af for 1 < k < 6 in Table 2.1. For more 
details, one can refer to [26]. 

3. Multistep schemes for decoupled 2FBSDEs. In this section, we first 
consider the multistep schemes for solving decoupled 2FBSDEs (2.9), i.e., the forward 
SDE is independent of (Y t , Z tl A t , T t ). To begin, let A be a positive integer. For the 
time interval [to,T], we introduce a regular time partition: 


to < ti < fjv — T. 

We set At tnj fe = t n+ k —t n for n, G {1, 2 ..., N} and k € N satisfying n + k < N. For 
t > t n , we denote A W tn ,k = W tn+k - W t „, A t tn , t = t-t n and AW t „, t = W t - W tn , 
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3.1. Four reference ODEs. Let © t = (X t . Y t , Z t ,A t ,F t ) be the solution of the 
decoupled FBSDEs (2.9) with terminal condition g(X?)■ We set E x ^ [•] = E [-| ^”’ x ] . 
By taking conditional expectation Ef^ [•] on both sides of the last two equations in 
(2.9), we obtain 

(3.1) El[Y t ]=El[g(X T )] +j\l[f(s,Q s )]ds, t £ [t n , T], 

(3.2) El [Z t ] =El [Z t J + jf El [A.] ds, t £ [t n , T], 

By taking derivative with respect to t in (3.1)-(3.2), one gets the following two refer¬ 
ence ODEs: 

(3.3) ^ = -EU/(*>e t) ], t&[t n ,T\, 

(3 - 4) =E i[A t]j te[t n ,T} 

Under sufficient regularity assumptions on the given data, the integrand E x n [/(s, 0 S )] 
and El [A t ] are continuous at s = t. Note that we also have 

Y tn =Y t + J^ f(s,@ s )ds - Z s dW s , t £ [t n ,T] 


By multiplying both sides of the above equation by A t (where ” T ” stands for 
the transposition), and taking the conditional expectation E x [■] on both sides of the 
derived equation, we obtain, for t £ [f„,T] 


(3.5) 0 = El [Y t AW t l t ] + J* El [f(s, 6 S )A^J ds - J* E^ [Z a ] 


ds, 


Similarly, for the last equation in (2.9), we have for t £ [t n ,T] 


(3.6) 0=E^ [ZJAWIj] 


El [A T s AW t T ntt ] ds 



[r s ] ds. 


Now, by taking derivative with respect to t £ [t n ,T) on both sides, one gets the 
following reference ODEs: 


(3.7) 

(3.8) 


dEf„ [Y t AWl tt ] 

dt 

dEl [ZJAWlj] 

dt 


~El [f(t,& t )AW t l t ] +E l[Z t ], 

El [.AJAW t l t ] +EI [r t ], 


t £ [ t n , T ] 
t £ [ t n , T] 


The equations (3.3), (3.4), (3.7), and (3.8) are reference ODEs for the decoupled 2FB- 
SDEs (2.9). Our multistep numerical schemes will be constructed by approximating 
the associate derivatives and the conditional expectations in these ODEs. 

3.2. The semi-discrete scheme. Motivated by Theorem 2, we choose smooth 
functions b(t,x ) and d(t,x) for t £ [t n ,T] and x £ R m with the constraints b(t n ,x) = 
b(t n ,x ) and d{t n ,x) = a(t n ,x). Let X* n,x be the diffusion process defined by 


(3.9) 


Xt 


t n ,X 


= x+ 6(s,X*”’ x )ds+ / d(s,Xl n ' x )dW t 
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Note that Y t ,Z t ,T t ,A t are all functions of (■ t,X t ). Now, let (Y t tn ' ,x 1 Zl n,x ) be the 
value of function ( Y tl Z t ) at the time-space point (t,Xl n ’ x ). , by Theorem 2, we have 


d£f„ [It] 

dEj n 

dE?„ [y t AlV- >4 ] 

dEl [*?»■* AW£, t ] 


dt 

t=t n dt 

dt 

t=t n 

dt 

t—t n 

t=t 

d£f„ [£t] 

d£? n 

dEl [Z?AW? nit ] 

dE L [(Zt^yAW t l 

,*] 

dt 

t=t B dt 

dt 

t—t n 

dt 

t=t n 


Now applying (2.12) to 


and 


dEfJ(Z^- x ) T AWT,,] 
d t 



dE?„[y t ‘— aw^.,1 

d^U^M 

dt 

ana dt 

’ dt 


t=t n 

t—t n 


, we deduce 

t=t n 


(3.10) 


Wt] 


d t 

dE^ [Y t AW t T n 



= 5> m e* 

t-t-n 2=0 
k 

t=t„ *=1 


t; 


+ 




I d« 


7t n ,X 


+ R k A,ni 


= E a ^ E " 

t—tn 2 = 0 

= J2a k ,El[(Zl:p T AWl 

t=t n *=i 


+ # 


r,n? 


Here are given by (2.11), and Ry tU , R^ n , R\ n and -^r n are the associate trun¬ 
cation errors, i.e. 


m = dE?„ [Y t ] 

xx y,n 


d t 


r>k 

Ix z,n 


TDK _ 

dt A n 


jDk _ 

n r,n ~ 


t—t n 2=0 

dEl [Y t AW t l, t ] 


-J2 a ^ E t 


y*n,z 

tn-\-i 


d t 

dE? [z t ] 


t=t n i=1 


Ytl X AWli 


d t 


t~t n 2=0 
dE?„ [Z t T A^ i( 


-E a ^ E " 


ytn,X 
^tn + i 


dt 


-J> M E* n [(Z^VAWY 

t=t~ *= i 
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By inserting (3.10) into (3.3), (3.4),(3.7) and (3.8), respectively, we deduce 


(3.11) 


k 


2 =0 

k 

Y tn+f ] = —f{t n ,x, Y tn ,Z tn ) + R 

2=1 

k 

Y tn+ AW^] 

— Zt„ + R k z ,n 

J2 a ^ E k 

2=0 

k 

An +i \ = At„ 

+ R h A,m 

J2 a ^ E k 

'zj n+i AWY 

II 

3*’ 

+ 

'-I5S- 

s 

2=1 




_ _ TDK TDK 

lx z,n'> Jri r,n 


_ TDk 


and R\ 


_ TDk 

— ~^Aa 


where R^ n = -R^ n ,R^ n 

Let Y n , Z n , A n and r n be four random variables that represent the approximate 
values of the solutions Y t , Z t , A t and T t of the 2FBSDE in (2.9) at time t n , respectively. 
By removing the truncation error terms Ry n , R\ n and R^ n from (3.11), we 

obtain the following semi-discrete numerical scheme for solving the 2FBSDEs (2.9): 


Scheme 1. Assume that random variables Y N 1 and Z N 1 , i = 0,1,..., k — 1 
are known. For n = N — k,... ,0, with X*”’ A being the solution of (3.9), solve 
Y n = Y n (X n ), Z n = Z n (X n ), A n =A n (X n ) audT n = T n {X n ) by 


k 


(3.12) 

z n = Y, a ^< 

3 =1 

k 

(3.13) 

r n = E^X 

3 =1 

k 

(3.14) 

3=0 

k 

(3.15) 

-ctk,o Y n = 5>,E a ; 
1 =i 


“ [Y n+ >XWY\ , 

“ [(Z n+ i) T AWY] , 

“ [Z n+j ],, 

“ [Y n+j ] + f(t n ,X n ,Y n ,Z n ,T n ), 


where Y n+ ^ and Z n+ i are the values of Y n+ i and Z n+ j at the space points X*"^ . 
In addition, we choose b and a in (3.9) as 

(3.16) b(s,Xl n ’ x ) =b(t n ,x), (r(s,X t s n ’ x ) = a(t n ,x), Vse[t„,T]. 

Then, (3.9) yields the Euler discretization scheme for the forward SDE. In this case, 
Scheme 1 becomes 

Scheme 2. Assume that random variables Y N ~ l and Z N ~ l , i = 0, 1,..., k — 1, 
are known. For n = N — k,... ,0, solve the random variables X n '° , (j = 1,2,..., k ), 
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Y n , Z n , 

A n and F n by 

(3.17) 

X n ' j =X n + b(t n ,X n )At n>j + cr(t n ,X n )AW n j, j 

(3.18) 

k 

Z n = P n+jAW n,i\ > 

3 = 1 

(3.19) 

A"=J2a k , i£" [Z^], 

3=0 

(3.20) 

F n = [(z n+j ) T A WW], 

3= 1 

(3.21) 

k 

-a k>0 Y n = [y n +i] + /(t T , ) x n ,y n ,z",r*), 

3=1 


where Y n+: > and Z n+ i stand for the values of Y n+ i and Z n+: > at the space points 
X n ’°, respectively. 

If the functions E* [Y t ],Ef n [Y t [Z t \, Ef n [Zf AW t T nA ] and their deriva¬ 

tives (with respect to t) up to order k + 1 are bounded, we have the estimates 

(3.22) Ry n = O (At) k , m n = 0(At)\ R k ^ n = O (At) k , R k , n = 0(At) k . 

3.3. The fully discrete scheme. We now propose our fully discrete schemes. 
To this end, we first introduce the time-space partition. Given a time partition T := 
{to, ti,..., tjv}, we introduce a series of space partitions T>h '■= {'E > h n }n=o,i,...,N for 
space R m , with each V^ corresponds to the time level i„eT and h = (ho, hi ,..., ftjv) 
being the partition density of each . We use V 7 ^ for simplicity if the context is 
clear. In addition, If the space partitions w.r.t. all the time levels are the same, then 
T>h is used to represent the unified space partition. 

For the space partition = {xj £ R m |j £ I, I C N}, Xj is called the grid 
points and / is called the index set of the grid. The partition density h n of is 
defined by h n = max d(xi,Xj ), with Xi,Xj adjacent in T> 7 f , where d(-,-) is the 

Xi.XnE'DV; n 

J n n 

distance between two points in M m . Given x £ R m , a neighbor gird U^ x of x, is a 
finite-element subset of VV: satisfying min d(x,Xi) < min d(x,Xi). 

hn *<evj! n VrZ. 

With the above spacial partition, we seek to solving Y n , Z n , A n and T n at grid 
point x £ T> r f . More precisely, for each x £ V^, n = N — k,..., 0, we aim at solving 
Y n , Z n , T" and A n by 

k 

3 =1 
k 

T n = ]To fe ,X„ [(Z n+ >) T AWX] , 

3 =1 
k 

An =J2 a ^ E i [ 2n+j ]’ 

3=0 
k 

-a k , 0 Y n ='}Ta k , j El [?"+'] + f(t n ,x,Y n ,Z n ,T n ), 

3 =1 


(3.23) 
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where Y n+ i and Z n+ ' J are respectively the values of Y n+ i and Z n+ i at the space 
points X n ’ ; i that is defined by 

(3.24) X n,j = x+ b(t n ,x)At n j + a(t n ,x)AW n j, j = l,...,k. 

Plug (3.24) into the conditional expectation of (3.23), since Y, Z,T and A are all 
functions of X nj , we have 

E?„ [Y n+j AWX] = E^ n [Y n ^(x + b(t n ,x)At n j+a(t n ,x)AW nJ )AWX] 

El [Z n+ q = El [Z tn+j (x + b(t n ,x)At nj +a(t n ,x)AW nj )\ 

El [(^ n+i ) T A^ T j =El [(Z n+ i) T (x + b(t n ,x)At nd +a(t n ,x)AW nJ )AWl j \ 
El [Y n +i] =El [ Y n+ i(x + b(t n ,x)At n j+a{t n ,x)AW nJ )] 

Note that quadrature methods should be applied when approximating the conditional 
expectations above. Any efficient quadrature rules such as the Monte-Carlo methods, 
the quasi-Monte-Carlo methods, and the Gaussian quadrature methods can be used. 
Here, we will use the Gaussian quadratures based on the zeros of Hermite polynomials, 
for details, one can refer to [26]. In what follows, we shall denote by E n ’ x [•] the 
numerical quadrature operator for the conditional expectations. 

Furthermore, when quadrature method is applied to approximate the conditional 
expectations, non-grid information might be needed. That’s to say, for x £ X>JJ, points 
(X nj defined by (3.24) e.t.c.) that are not in 'D r }', +J are needed when approximate 
the conditional expectation at time level t n +j. Thus, interpolation methods are also 
needed. We shall denote by Id, a local interpolation operator, such that lug is the 
continuous function interpolated by the values of g at the grid points V. We also 
denote by Ip with the superscript n indicating the interpolation method at time level 
t n . If the grid is a neighbor grid U^ x of x, we denote by I^, x instead of 1 ^, , for 
simplicity. Note that any interpolation methods can be used here, however, care 
should be made if one wants to guarantee the stability and accuracy. 

Now by introducing the operators E n,x [•] and K, (or more precisely I" + it„,* ) 

t n+j 

one can rewrite (3.11) in the following equivalent form 


k 


Z tn 

= E“W ^ 

3 =1 

^% n ,*Y tn+j AWX 

’ *n+j 

_ jxk ofc,E jxk ,I 

-^z^n ' lx z,n ' rL z,n’> 


k 





= ^> fej E^ 

3 =1 

L *n+j 

- K V n + H T n + H Vn 


k 





= 5> fc , 

3=0 

l v + k-' z *-« 

p/c i pk,E , jjk, I 

-*^4,71 ' A,n ' 


k 

p , 



&kfiYt n 

= ]Ta feJ E^ 

3 =1 

I n +3 y- 

L ’ t n+j J 

+ f(t n , x, Y tn ,Z tn ) 


Tjk 1 j^kj 
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where 


d/c, E 

-^2,71 

d/c, I 

1L~ y . 


R 


kj 


= Ej-1 - E".-) awv], 

= n; + j. 


,fc,E 

'A,n 

t fc,I 

A,n 

k, E 


-£) j£ z n > x ^ ti+j 
’ *n+i 


) A ^- 


= E‘=i«fc,i(E? B 7 ET ‘’ X ) [^-J» 


— Ej=l 

■\k 


x,h 


^tn+j ) 




i?; 

pfc,E 

-ft,,’ 


= EJ=i«fc,i(Ef n -E” lX ) 

= E- =1 « fej Ee [(^-1^3 
= - ELi - E n> ") [?" +j ], 


^ =-eL^X; 


7 T 

T>,Xj n ’ X tn+j 

n+j 


) A ^ 


Y t . 


— ll ^l +- , V 
u 75 X tn ' x 1 tn +i 

’ *"+3 


The four terms and i?p® are the error terms resulted from approxi¬ 

mating conditional expectations, and the other four terms Ry’ l nl ,, and i? r ’^ 
are the error terms caused by numerical interpolations. 


By removing those error terms Ry n , R k i?: 




,fc,i 




1? 


fc,E 


fc,I 

Xj y,n'> 


TDk TDk, E E>/c,I Di 

-* n • tl r n > 1 L ~ n i ' I 


2,n7 1L z,m 1L z,m lh A,n’> 1L A,ni 


R 


k, E 


and i? r ’ ra from (3.25), we propose our fully discrete scheme for 


A,n’ *T,n’ *T,n 

solving 2FBSDEs as follows: 

Scheme 3. Assume random variables Y N ~ l and Z N ~ l defined on D?~ l , i = 
0,1,..., A: — 1, are known. For n = N — k,..., 0, and for each x € 2?^, solve X n , 7", 
Z n , A” and T n by 


(3.26) 

X n ’ 3 = x + b(t n ,x)At I 

k 

(3.27) 

z n = J2 ak ’$‘ 

3= 1 

k 

(3.28) 

r n = E“^ 

3=1 

k 

(3.29) 

&. 

ii 

e 

(3.30) 

-a kfi Y n = 

j =i 


7 T ■ 


Tf^+l yn+j 
T5,X"J 


+ /(t ni s,y n ,z n ) r"). 


In Scheme 3, for a fixed integer k, to get the values Y n , Z n , A n and T n on each 
time level t n at the grid V^, we shall solve the values at x € T>% one by one. Since 
they are independent with each other, the precess can be completely parallel. For 
a fixed x G T>^, in Scheme 3 , we firstly solve X n ^ by the Euler scheme (3.26) for 
1 < j < k: and then solve Z n , T n and A n by (3.27), (3.28) and (3.29) explicitly; 
finally, solve Y n by (3.30) implicitly. Thus some iteration methods are required for 
solving Y n . 
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If the function f(t n , x , y , z, 7 ) is Lipschitz continuous with respect to y , for small 
time partition step size Af n , we can use the following iteration procedure to solve Y n 


(3.31) a M r l+1 = -^ aiJ r 

1=1 


f n +l yn+J 




with the iteration stopping condition \Y n ’ l+1 — Y n ' l \ < eo, where eo > 0 is a given 
tolerance. In case the function f(t n ,x,y,z, 7 ) is differentiable with respect to y , to 
accelerate the convergence rate, the Newton iteration method can be applied , i.e. 


a kfi Y n ’ 1 + £ a kJ E' 

1 _ Y n il __ 


]r n +J 


+ /(t„,a;,y n > , ,z",r") 


afe,0 + fy{tn, X , Y n ’ l ,Z n , r n ) 


Note that the local truncation errors of Scheme 3 consist of twelve terms „, 

p/c,E jDk, I pfc d/c,E p/c,I jDk p&,I d/c j rpi r 

ri y,n’> Ix y,n^ Ix z,n’) Ix z,n 5 1X1 z,n 5 rx A,n’> ri A,n'> rx A,n'> Jrx r,ni i *T,n cllia - r *T,n‘ x lie iUUi 

terms Ry^ n , R k . n: Ra,u an d ^r,« defined respectively in (3.10) come from the ap¬ 
proximations of the derivative, and the four terms R k f n , R k 'li , and defined 

in (3.25) are the local interpolation errors. For these eight terms, when data b , a, f 
and g are smooth, the following estimates hold (provided that degree r polynomials 
interpolation is used) 

(3.32) R k ytn , R% n , R% n , R k >n = 0(At) k ), R k / n , R k ’* n , R k / n , R*f n = O ( h r+l ) . 

The other four terms R~’ni R^An and Rr ^ are the local truncation errors re¬ 

sulted from the approximations of the conditional expectations. Noticed that these 
conditional expectations are functions of Gaussian random variables, thus, one can 
construct efficient Gauss-Hermite quadrature rule for their approximations. 

4. Extensions to coupled 2FBSDEs. In this section we extend Scheme 3 to 
solve the fully coupled 2FBSDE. More precisely, we first update Scheme 3 into the 
following scheme 

Scheme 4. Assume random variables Y N ~ l and Z N ~ l defined on ~D^ 1 , i = 
0,1,..., ft — l, are known. For n = N — k ,..., 0, and for each x £ Vsolve Y n , Z n , 
A n and T" by 

X n 'i =x + b(t n , x, Y n , Z n , r n )At nJ + a(t n , x, Y n ,Z n , T n )AW n>J , 

j = 1,2,..., ft, 

fc 

V'AI r,,^ 
i=i 
k 


(4.1) 


r n = E a wE" 

3 =1 
k 

A " = 5><bl E" 

3=0 

k 

-a k , 0 Y n =J2<*k,jE n 


K + LAz n+j V&w; 


]J n +l gn+j 


T),X n ’0 


Tf n +1 yn+j 
V,X n ’i 


+ f(t n ,x,Y n ,Z n ,T n ). 


1=1 
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The main difference in Scheme 4 is that X n, I, Y n ,Z n and T™ are all coupled 
together. Thus, one needs to solve nonlinear equations with some iterative procedure. 
In our numerical experiments, we use the following iterative scheme. 

Scheme 5. Assume random variables Y N ~ l , Z N ~ l and T N ~ l defined on , 
i = 0,1,..., k — 1, are known. I is the iterative time, Y n ’, Z n ’ 1 , A 71 ’ 1 and T rl ’ i are the 
corresponding values at the l-th iterative. For n = N — k,, 0, and each x G T) 7 ^, 
solve Y n , Z n , A n and T n by 

1. Initialize Y n ’°, Z n ’° and T 11 ’ 0 by Y n+1 , Z n+1 and T n+1 . separately, i.e. 

Y n ’° = Y n+1 (x), Z n ’° = Z n+1 (x ) and T n ’° = T n+1 (a;); 

2. Calculate Y n ’ l+1 , Z n ’ l+1 , A n ’ l+1 and T n ’ i+1 fori = 0,1,... by 


X n ’i =x + b(t n , x, Y n ' 1 , Z n ’ 1 , T"> , )At nJ + a{t n , x, Y n ’ 1 , Z n \ T n ’ l )AW ntj , 
j = 1 , 2 ,.. ,,k, 

k 

Z n,i+1 = J2a kj t n ’ x t + ] n / n+J AH/h 


l=i 

k 


r n,;+i = E n ’ x l£i nii (Z n+j ) T AWY 

3 =1 
k 


H n +3_ yn+j 
V,X n >3 


A n,i+i = J2 a kJ E n 
l=o 
k 

-a k , 0 Y n ’ l+1 =^a k jE n ’ x |l” + j-„, 3 -F n+J I + f(t n , x, Y n ’ l+1 , Z n ’ l+1 , T"’ i+1 ) 


l=i 


Repeat the above calculation, until 

max { \Y n ’ l+1 - Y n ’ 1 1, | Z n ' l+1 - Z n ’ l \, \A n ’ I+1 - A n ' l \ , |T n ’ i+1 - r n ’ ; |} < e 0 ; 

3. The values (Y n , Z n ,T n , A n ) are set to be (Y n ’i +1 , A n ’ l+1 , T"’ i+1 ). 

The computation is complete. 

It is worth to remark that our Schemes 4 and 5 are heuristic generalizations of 
the decoupled ones in the last section. In case the drift coefficient b and the diffusion 
coefficient er do not depend on Y, Z and T, Scheme 5 coincides with Scheme 3. 

Remark 4.1. We have finished the construction of our multistep numerical 
schemes , i.e., Scheme 3 for decoupled 2FBSDEs and Scheme 5 for coupled 2FBSDEs. 
In both cases, the Euler method is used to discrete the forward SDEs, thus, the total 
computational complexity can be significantly reduced. Furthermore, we shall show by 
numerical tests that the solution {Y t , Z t ,T tl A t ) can still be of high order accuracy, in 
the next section. 


5. Numerical experiments. In this section, we shall provide with several con¬ 
structive numerical tests to show the efficiency and accuracy of our multistep schemes 
proposed in the previous sections. Applications of our numerical methods for stochas¬ 
tic optimal control problems will also be presented. 

In all our numerical tests, we shall consider the uniform partition for the time- 
space domain [0,T] x W n , for simplicity. That is, for a given TV > 0,h > 0, we 
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set 

T := {t n | t n = nAt, n = 0,1,. .., TV, At = ^}, 

T> h ■= {*j \ Xj = x 0 +i ■ h, j := ■ ■ ■ , jm) T , for each j t £ Z}, 

with xq £ K. m the initial state. Specially, for the one dimensional case, i.e. m = 1, 
this yields T>h = {xj \ Xj = xq + jh,j = 0, ±1, ±2,...}. Our scheme is applied to 
calculate the values of Y n , Z n , T” and A n , which are the numerical values of Y, Z, 
r and A at every time-space point (t n ,x$) £ T x T>h■ 

In our tests, we shall adopt the local Lagrange interpolation method Id, so that 
the interpolation error estimates in (3.32) holds, and the Hernrite-Gaussian quadrature 
rule will be used to approximate the associated conditional expectations. Since we 
aim to checking the convergence rates of our multistep schemes, all errors resulted 
from these approximations are controlled by choosing appropriate parameters (e.g. 
interpolation order r and the numerical of Gaussian-Hermite points). In particular, 
we shall use 10 Gauss-Hermite quadrature points in each dimension, and for the 
Lagrange interpolation, more than 6 points will be used. To balance the errors result 

fc + 1 

from the time discrete truncation and the space truncation, we choose h = (A t) r+1 , 
where r is the degree of the Langrange interpolation polynomial. Furthermore, for 
the fc-step schemes, the information for {Y N ~i, Z N ~^}j =1 are needed, this can be 
obtained by using other numerical methods with small time steps, to maintain the 
high order convergence rates. Here, we just artificially assume that the values are 
known for simplicity. 

The numerical results are obtained with FORTRAN 95 on a workstation with one 
Intel Xeon E5-2620 v2 CPU (12 cores, 2.10 GHz ). To accelerate the perfor¬ 
mance, OpenMP techniques are used. To guarantee the computing precision, we use 
the long double type (real (16)) digital for the float variables when programming. 
The long double variable has 34 significant digits which supplies enough precision for 
our computation. However, the cost of time is increasing dramatically compared with 
variables all defined as double (real (8)), but even so, the time elapsed of our pro¬ 
grams is still competitive. In what follows, we will denote by CR the convergence rates 
and T r the running time, respectively. For all our numerical examples, the terminal 
time T is set to be 1.0. 

5.1. Decoupled and coupled 2FBSDEs. We first test Scheme 3 for solving 
decoupled 2FBSDEs. The first example considered is 

{ dX t = sin (t + X t )d t + ccos(f + X t )d W tl 

-d Y t = ( - cos(f + X t )-Z t - cos(f + X t )(Y t 2 + Y t ) - ir t )dt - Z t dW t , 

dZ t = A t dt + r t d W t , 

with the initial condition X to = x and the terminal condition Yp = sin(T + Xp)- It 
can be shown that the exact solutions are 

Y t =sm(t + X t ), Z t = ccos 2 (f + X t ), T t = -2c 2 sin(i + X t ) cos 2 (t + X t ), 

A t = — csin(2f + 2X t )(l + sin(t + X t )) — c 3 cos(2 1 + 2X t ) cos 2 {t + X t ). 

In the numerical test, we set x = 0.5 and to = 0, and solve the 2FBSDEs by scheme 3 
with different step parameter k. The numerical errors IF 0 — To I , I Z° — Zq I , IF 0 — Fo I 
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and | A 0 — A 0 | and the corresponding convergence rates are listed in Table 5.1. It is 
shown in Table 5.1 that: the multistep numerical scheme works very well, and the 
numerical error goes to machine accuracy. Moreover, the method admits a /c-order 
convergence rate, and it remains stable for 1 < k < 6, which is coincide with the 
classic numerical ODEs theory and our previous results [26]. 

Furthermore, given a fixed accuracy tolerance, it would be more efficient if the 
multistep scheme with a large k is used. This can be more easily seen from the 
following table. 


STEP 

N 

T rn 

O 

2 U -Z 0 

| -| 

o 

n 

o 

O 

1 

o 

k = 1 

8192 

135.0s 

5.172E-05 

2.739E-05 

7.573E-07 

3.585E-07 

CM 

II 

-se 

2048 

10.82s 

1.232E-09 

4.704E-08 

2.301E-08 

1.920E-07 

CO 

II 

512 

5.58s 

3.423E-09 

3.097E-09 

7.797E-10 

4.677E-09 

k = 4 

128 

2.85s 

1.159E-09 

1.733E-09 

2.113E-09 

6.042E-09 


We now consider an example with the geometry Brownian motion in the forward 
SDE, the main feature here is that the drift and diffusion terms are unbounded. The 
example yields 


( dX t = rX t dt + cX t dW u 

( 5 . 2 ) | -d Y t = -e-% + y^ 2 Y t -y + ^Z t dt - Z t dW t , 

l dZt = Afdt + T t dW t 


with terminal condition Yt = Te . It can be shown that the exact solution is 

4c 2 


Y t = te~ 


z, = 

M 4 


A t = ~Jf e " K 1 + 2rt + c2 ^ X t ~ tX t 


r 2„ — 

2 c 2 t 
~M~ 


2 r + c 2 




In this example, we set Xo = 1.5, and the parameters are chosen as r = 0.2, c = 
0.01, M = 4. The corresponding numerical results are listed the in Table 5.2. Again, 
our multistep schemes behaves very well, and high order convergence rates are ob¬ 
tained. 

We now test Scheme 5 for solving coupled 2FBSDEs. Here, iterative process is 
needed. To show the high accuracy of the proposed schemes, we set the tolerance as 
e = 10“ 25 . We first consider the following coupled 2FBSDEs: 


dX t = (sin(f + X t ) + Z t /c + sin(t + X t )Y t - l)dt 

+ (ccos (t + X t ) - c + ccos 2 (f + X t ) + csin(f + X t )Y t )dW t , 

— d Y t = — cos(t + Xt)—Zt — cos {t + X t ){Y t 2 + Y)) — —Tjdt — ZtdWt, 
c 4 

, dZt = Atdt + r t dWf, 


with conditions Yt = sin(T + Xt) and X to = x. The exact solution is 

Y t = sin(i + X t ), Z t = ccos 2 (t + X t ), r t = — 2c 2 sin(i + X t ) cos 2 (t + X t ), 
A t = —csin(2t + 2X t )(l + sin(t + X t )) — c 3 cos(2 1 + 2X t ) cos 2 {t + X t ). 
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Table 5.1 

Numerical results of example 1 with c = 0.1 


STEP 


N 

C 

1 

5^ 

|2 u -Zo| 

"I 

C 

1 

•n 

o 

U- 

C 

1 

o 

Trn 

K=1 

DT 

32 

1.253E-02 

6.800E-03 

3.156E-04 

1.703E-03 

1.87s 

64 

6.411E-03 

3.466E-03 

1.854E-04 

1.059E-03 

128 

3.248E-03 

1.752E-03 

8.667E-05 

4.409E-04 

256 

1.637E-03 

8.808E-04 

3.570E-05 

1.381E-04 

512 

8.227E-04 

4.405E-04 

1.411E-05 

3.472E-05 

CR 


0.98 

0.99 

1.13 

1.42 

K=2 

DT 

32 

6.123E-06 

1.897E-04 

1.086E-04 

9.312E-04 

4.96s 

64 

1.334E-06 

4.761E-05 

2.614E-05 

2.210E-04 

128 

3.332E-07 

1.192E-05 

6.207E-06 

5.196E-05 

256 

8.494E-08 

2.995E-06 

1.500E-06 

1.253E-05 

512 

2.133E-08 

7.514E-07 

3.707E-07 

3.098E-06 

CR 


2.03 

2.00 

2.05 

2.06 

K=3 

DT 

32 

1.149E-05 

1.409E-05 

2.846E-06 

1.988E-05 

13.96s 

64 

1.615E-06 

1.668E-06 

3.734E-07 

2.382E-06 

128 

2.117E-07 

2.026E-07 

4.848E-08 

2.972E-07 

256 

2.704E-08 

2.496E-08 

6.202E-09 

3.727E-08 

512 

3.415E-09 

3.096E-09 

7.818E-10 

4.650E-09 

CR 


2.93 

3.04 

2.96 

3.01 

K=4 

DT 

32 

2.853E-07 

4.169E-07 

6.203E-07 

2.171E-06 

17.74s 

64 

1.843E-08 

2.700E-08 

3.588E-08 

1.122E-07 

128 

1.161E-09 

1.736E-09 

2.097E-09 

5.886E-09 

256 

7.271E-11 

1.106E-10 

1.275E-10 

3.463E-10 

512 

4.545E-12 

6.973E-12 

7.852E-12 

2.097E-11 

CR 


3.99 

3.97 

4.07 

4.17 

K=5 

DT 

32 

1.124E-08 

3.151E-08 

2.563E-08 

1.660E-07 

34.41s 

64 

5.297E-10 

9.433E-10 

8.493E-10 

4.027E-09 

128 

1.877E-11 

2.770E-11 

2.632E-11 

9.995E-11 

256 

6.180E-13 

8.174E-13 

8.310E-13 

2.951E-12 

512 

1.971E-14 

2.480E-14 

2.693E-14 

9.838E-14 

CR 


4.80 

5.07 

4.97 

5.18 

K=6 

DT 

32 

4.721E-10 

9.254E-10 

5.225E-09 

1.144E-08 

50.75s 

64 

7.847E-12 

1.542E-11 

7.602E-11 

1.545E-10 

128 

1.144E-13 

2.506E-13 

1.140E-12 

2.144E-12 

256 

1.741E-15 

4.074E-15 

1.710E-14 

2.863E-14 

512 

2.674E-17 

6.557E-17 

2.622E-16 

4.232E-16 

CR 


6.03 

5.94 

6.06 

6.18 

K=7 

DT 

32 

1.422E-11 

1.805E-11 

4.311E-10 

7.197E-09 

83.65s 

64 

3.876E-13 

8.371E-13 

1.417E-12 

3.941E-12 

128 

2.365E-15 

1.959E-15 

2.089E-14 

3.823E-13 

256 

2.762E-17 

1.366E-17 

1.265E-15 

1.528E-14 

512 

4.993E-18 

1.958E-16 

9.029E-17 

1.894E-16 

CR 


5.67 

4.89 

5.45 

5.84 

K=8 

DT 

32 

7.015E-12 

1.501E-10 

2.182E-09 

2.147E-08 

108.42s 

64 

4.860E-12 

6.270E-11 

3.480E-10 

7.047E-09 

128 

1.117E-10 

6.105E-09 

8.660E-08 

1.014E-06 

256 

7.367E-03 

3.798E+00 

4.029E+02 

3.354E+04 

512 

NaN 

NaN 

NaN 

NaN 

CR 


NaN 

NaN 

NaN 

NaN 


The numerical results are shown in Table 5.3, where one can concludes that the 
numerical methods work well and high-order convergence rates are obtained. However, 
due to the computational complexity, numerical results with only multistep up to 
k = 3 are presented. 
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Our next couple 2FBSDEs example yields 

1 


Xt — x + 
Y t = 


q (1 + exp (s + X s ))(l + Y s ) J q 

exp (T + Xt) 


ds+ / V s dlF s , 


1 + exp (T + X T ) 


2 Y„ 


1 


r s - 


y,z s 


1 + 2 exp (s + X s ) 2 \ 1 + exp (s + X s ) J 

z t = z T - [ A s ds- / r s div s . 


ds — 


Z s dW s , 


Table 5.2 

Numerical results of example 2. 


STEP 


N 

C 

1 

\Z U — Zo\ 

O 

u 

1 

3 

C 

1 

X- 

o 

Tm 

K=1 

DT 

32 

4.808E-03 

2.331E-04 

2.431E-06 

6.699E-07 

2.58s 

64 

2.402E-03 

1.169E-04 

1.232E-06 

3.335E-07 

128 

1.200E-03 

5.854E-05 

6.207E-07 

1.665E-07 

256 

5.999E-04 

2.929E-05 

3.115E-07 

8.308E-08 

512 

2.999E-04 

1.465E-05 

1.560E-07 

4.143E-08 

CR 


1.00 

1.00 

0.99 

1.00 

K=2 

DT 

32 

1.385E-05 

1.294E-06 

1.459E-07 

1.457E-06 

7.54s 

64 

3.582E-06 

3.252E-07 

3.653E-08 

3.603E-07 

128 

9.098E-07 

8.147E-08 

9.140E-09 

8.954E-08 

256 

2.291E-07 

2.039E-08 

2.286E-09 

2.231E-08 

512 

5.747E-08 

5.098E-09 

5.717E-10 

5.569E-09 

CR 


1.98 

2.00 

2.00 

2.01 

K=3 

DT 

32 

3.490E-07 

5.124E-08 

1.213E-09 

8.537E-09 

16.73s 

64 

4.638E-08 

6.484E-09 

1.552E-10 

1.130E-09 

128 

5.983E-09 

8.157E-10 

1.960E-11 

1.454E-10 

256 

7.604E-10 

1.023E-10 

2.461E-12 

1.846E-11 

512 

9.590E-11 

1.281E-11 

3.082E-13 

2.326E-12 

CR 


2.96 

2.99 

2.99 

2.96 

K=4 

DT 

32 

4.889E-09 

8.085E-10 

7.806E-11 

6.895E-10 

34.39s 

64 

3.369E-10 

5.171E-11 

4.934E-12 

4.232E-11 

128 

2.212E-11 

3.276E-12 

3.100E-13 

2.601E-12 

256 

1.417E-12 

2.063E-13 

1.942E-14 

1.610E-13 

512 

8.969E-14 

1.295E-14 

1.215E-15 

1.003E-14 

CR 


3.94 

3.98 

3.99 

4.02 

K=5 

DT 

32 

4.916E-11 

7.997E-12 

6.357E-14 

2.573E-12 

61.72s 

64 

1.548E-12 

2.514E-13 

1.625E-15 

7.997E-14 

128 

4.545E-14 

7.896E-15 

3.486E-17 

2.277E-15 

256 

1.316E-15 

2.481E-16 

6.573E-19 

6.360E-17 

512 

3.843E-17 

7.792E-18 

1.067E-20 

1.790E-18 

CR 


5.08 

4.99 

5.63 

5.12 

K=6 

DT 

32 

2.246E-12 

3.449E-13 

3.677E-14 

4.326E-13 

93.42s 

64 

3.805E-14 

5.534E-15 

5.708E-16 

6.289E-15 

128 

5.870E-16 

8.807E-17 

9.179E-18 

9.283E-17 

256 

8.808E-18 

1.394E-18 

1.463E-19 

1.391E-18 

512 

1.319E-19 

2.197E-20 

2.318E-21 

2.136E-20 

CR 


6.01 

5.98 

5.98 

6.07 

K=7 

DT 

32 

3.908E-13 

1.027E-13 

2.830E-15 

1.527E-13 

196.54s 

64 

4.974E-16 

1.217E-15 

3.787E-17 

1.929E-15 

128 

3.750E-17 

1.204E-17 

3.601E-18 

1.489E-16 

256 

3.893E-19 

2.617E-17 

2.460E-15 

3.265E-13 

512 

5.287E-17 

8.003E-15 

6.419E-12 

1.677E-09 

CR 


3.60 

1.29 

-2.83 

-3.42 




























































18 


T. Kong, W. Zhao, and T. Zhou 


Table 5.3 

Numerical results of example 3. 


STEP 


N 

o 

1 

$ 

\Z U - Z 0 \ 

O 

Eh 

1 

3 

u_ 

c 

1 

o 

Tm 

K=1 

DT 

32 

4.126E-02 

9.266E-03 

3.980E-03 

5.001E-02 

35.42s 

64 

2.085E-02 

4.782E-03 

1.996E-03 

2.532E-02 

128 

1.048E-02 

2.430E-03 

9.968E-04 

1.273E-02 

256 

5.255E-03 

1.225E-03 

4.977E-04 

6.377E-03 

512 

2.631E-03 

6.149E-04 

2.486E-04 

3.192E-03 

CR 


0.99 

0.98 

1.00 

0.99 

K=2 

DT 

32 

5.077E-04 

4.045E-05 

2.561E-04 

1.640E-03 

232.23s 

64 

1.250E-04 

7.443E-06 

6.648E-05 

4.266E-04 

128 

3.099E-05 

1.541E-06 

1.690E-05 

1.086E-04 

256 

7.713E-06 

3.462E-07 

4.258E-06 

2.740E-05 

512 

1.924E-06 

8.166E-08 

1.069E-06 

6.876E-06 

CR 


2.01 

2.23 

1.98 

1.98 

K=3 

DT 

32 

6.995E-05 

2.551E-05 

2.388E-05 

1.568E-04 

987s 

64 

8.908E-06 

3.156E-06 

2.693E-06 

1.840E-05 

128 

1.123E-06 

3.923E-07 

3.173E-07 

2.225E-06 

256 

1.409E-07 

4.889E-08 

3.843E-08 

2.733E-07 

512 

1.764E-08 

6.103E-09 

4.726E-09 

3.383E-08 

CR 


2.99 

3.01 

3.07 

3.04 


The associated exact solution is 

_ exp (t + X t ) _ exp (t + X t f 

1 l + exp(f + X t )’ t (1 + exp (t + X t )) 3 ’ 

_ exp(t + X t f(2 -exp(t + X t )) 

(1 + exp (t + X t )) 5 

A _ 2 exp (t + X t f (2 - exp (t + X t )) 

1 (1 + exp (t + X t )) 3 (l + 2 exp (t + X t )) 

exp (t + X t ) 4 (exp (t + X t ) 2 - 7exp (t + X t ) + 4) 

2(1 + exp (t + X t )) 7 

The corresponding numerical results are listed in Table 5.4. Similar convergence 
results are shown as the above numerical tests. 

5.2. Applications to stochastic optimal control. We now show that one can 
solve stochastic optimal control problems in a 2FBSDEs way by using our multistep 
numerical schemes. To this end, let us consider the control problem, whose dynamic 
state equation is described by a forward SDE 

(5.4) dX t = b(t, X t , a t )dt + cr(t, X t , a t )dW t , 
with the cost functional 

(5.5) J{a)=E 

The goal is to minimize the cost functional, i.e. Find a* £ U (where U contains all 
admissible controls, see e.g. [16] for the corresponding definition) satisfying 

V(t,x) := J{t,x\a*) = inf J(t,x;a). 

cx.e.U 


f(t, Xt,at)dt + g(X t ) 
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Table 5.4 

Numerical results of example 


STEP 


N 

o 

1 

$ 

\Z U - Z 0 \ 

O 

Eh 

1 

3 

u_ 

C 

1 

X- 

o 

Tm 

K=1 

DT 

32 

1.079E-03 

2.877E-03 

3.847E-03 

3.427E-03 

21.54s 

64 

5.374E-04 

1.463E-03 

2.008E-03 

1.796E-03 

128 

2.680E-04 

7.389E-04 

1.030E-03 

9.216E-04 

256 

1.338E-04 

3.715E-04 

5.225E-04 

4.678E-04 

512 

6.681E-05 

1.864E-04 

2.634E-04 

2.359E-04 

CR 


1.00 

0.99 

0.97 

0.97 

K=2 

DT 

32 

5.131E-05 

1.750E-05 

2.281E-04 

3.383E-04 

70.10s 

64 

1.328E-05 

4.871E-06 

5.814E-05 

8.806E-05 

128 

3.388E-06 

1.248E-06 

1.477E-05 

2.276E-05 

256 

8.523E-07 

3.171E-07 

3.691E-06 

5.698E-06 

512 

2.139E-07 

7.953E-08 

9.282E-07 

1.432E-06 

CR 


1.98 

1.95 

1.99 

1.97 

K=3 

DT 

32 

3.134E-06 

6.034E-06 

4.653E-06 

2.167E-05 

259.11s 

64 

4.058E-07 

7.806E-07 

5.391E-07 

2.879E-06 

128 

5.121E-08 

9.992E-08 

5.827E-08 

3.605E-07 

256 

6.407E-09 

1.257E-08 

7.059E-09 

4.530E-08 

512 

8.052E-10 

1.574E-09 

8.454E-10 

5.398E-09 

CR 


2.98 

2.98 

3.11 

2.99 


One can construct the corresponding HJB equation as 

f) C rr(t t cv') 2 c) 2 D 1 

^V(t, x ) + inf | 2 ^V(t, x) + b(t, x, x) - f(t, x,a)j= 0 . 

and moreover, we have 

*n \ ■ r f <x(t,x,a) 2 d 2 V dV \ 

(5.6) a (t, x) = arg mf j + b(t, x, a) — - f(t, x,a)\ . 

By inserting (5.6) into the BJB equaion, one shows that the cost function satisfies 

(5.7) ^V{t,x) +G(t,x, -^-V(t,x), V(t,x)) = 0, 

where 

(j(f 'T* ry* 

G{t, x,p, P) = - L y-- P + b(t, x, a*)p - f(t, x, a*). 

For the above nonlinear PDEs, one can construct a corresponding 2FBSDEs. 
Then, one can solve the associate 2FBSDEs using our numerical schemes to obtain 
(Xt, Y t , Z t , r t ), which yields the optimal state for the control problem, and finally, we 
can get obtain in view of (5.6) that 

(5.8) a* t =g(X t ,Y t ,Z t ,T t ) 

with g(-) being certain functions. This is a new approach dealing with the optimal 
control problem in a 2FBSDEs way. Now, we illustrate the idea by the following 
example. 

Tracking a particle under the microscope: consider the following system 

(5.9) dX t = (3atdt + adWt, 

where X t is the distance between the particle and the focus of the microscope, /3 € R 
is the gain in our servo loop and a > 0 is the diffusion constant of the particle. We 
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would like to keep the particle in focus, i.e. we expect that X t is as close to zero as 
possible. However, we have to introduce a power constraint on the control as well, 
as we cannot drive the servo motor with arbitrarily large input powers. We thus 
introduce the control cost 


J(a) = E 


p f X 2 dt + q f a 2 At 
Jo Jo 


where p, q > 0 allows us to select the tradeoff between good tracking and low feedback 
power. One can construct the associate bellman equation 

0 = + fc^V(t,x) + i>* 2 + „a 2 } 

= - ^{Xv(t,x)?+ W > 

with V (T, x) = 0. Furthermore, one can show that the optimal control parameter 
(i.e., exact solution) is 

(5.!°) ai t = -]LJLv(t,x). 

The classic numerical solution would relies on solving the above bellman equation. 

We now solve the problem in a 2FBSDEs way, to this end, we first construct the 
corresponding 2FBSDEs as follows 

{ AX t = (3cdt + crdW t , 

—d Y t = (~^Z 2 - —Z t + P X 2 )dt - ZtdW t , 

4 qa z a 

d Z t = A t dt + T t dW t . 


By solving this 2FBSDEs, we obtain the numerical solution in view of (5.10) by 

(5.12) a n = —Z n . 

2 qa 

In the numerical test, we set p = 0.1, r = 0.03, cr = 0.5, c = 0.1. The numerical 
results are shown in Table 5.5. It can be seen from Table 5.5 that the approach is of 
high order accuracy, both for the 2FBSDE solution and the optimal control a. 

6. Conclusions. We have extended our multistep schemes in [26] to the multi- 
step schemes 3 and 5 for solving the second order FBSDEs. The key feature of the 
proposed multistep schemes is that the Euler method is used to discrete the forward 
SDE, which dramatically reduces the entire computational complexity. Furthermore, 
it is shown that the quantities of interest (e.g., the solution tuple (Y t , Z tl A t ,T t ) in 
the 2FBSDEs) are still of high order accuracy. Several numerical examples are pre¬ 
sented to show the effective of the proposed numerical schemes. Applications of our 
numerical schemes for stochastic optimal control problems are also discussed. 

There are, however, some other related topics that need to be investigated: 

• High dimensional problems. Note the methods here can be easily extended to 
high dimensional problems. However, we have proposed the local Lagrange 
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Table 5.5 

Numerical results of stochastic optimal control. 


STEP 


N 

o 

1 

£ 

\Z U -Z 0 \ 

>1 

C 

1 

"I 

o 

|o! u — a* | 

T vn 

K=1 

DT 

32 

5.647E-03 

5.856E-02 

1.420E-02 

1.952E-02 

2.05s 

64 

2.788E-03 

2.912E-02 

7.075E-03 

9.706E-03 

128 

1.385E-03 

1.452E-02 

3.531E-03 

4.839E-03 

256 

6.903E-04 

7.249E-03 

1.764E-03 

2.416E-03 

512 

3.446E-04 

3.622E-03 

8.814E-04 

1.207E-03 

CR. 


1.01 

1.00 

1.00 

1.00 

K=2 

DT 

32 

1.438E-03 

8.777E-04 

5.823E-05 

2.926E-04 

6.74s 

64 

3.660E-04 

2.192E-04 

1.471E-05 

7.308E-05 

128 

9.230E-05 

5.478E-05 

3.696E-06 

1.826E-05 

256 

2.317E-05 

1.369E-05 

9.264E-07 

4.564E-06 

512 

5.806E-06 

3.422E-06 

2.319E-07 

1.141E-06 

CR 


1.99 

2.00 

1.99 

2.00 

K=3 

DT 

32 

3.439E-06 

4.495E-06 

1.848E-06 

1.498E-06 

19.44s 

64 

4.758E-07 

5.315E-07 

2.276E-07 

1.772E-07 

128 

6.243E-08 

6.456E-08 

2.824E-08 

2.152E-08 

256 

7.992E-09 

7.953E-09 

3.517E-09 

2.651E-09 

512 

1.011E-09 

9.869E-10 

4.387E-10 

3.290E-10 

CR 


2.94 

3.04 

3.01 

3.04 

K=4 

DT 

32 

4.583E-07 

3.240E-07 

3.570E-08 

1.080E-07 

57.62s 

64 

2.986E-08 

2.005E-08 

2.264E-09 

6.682E-09 

128 

1.902E-09 

1.246E-09 

1.425E-10 

4.154E-10 

256 

1.200E-10 

7.767E-11 

8.933E-12 

2.589E-11 

512 

7.535E-12 

4.847E-12 

5.592E-13 

1.616E-12 

CR 


3.97 

4.01 

3.99 

4.01 

K=5 

DT 

32 

4.302E-09 

2.299E-10 

8.347E-10 

7.662E-11 

94.61s 

64 

1.512E-10 

1.904E-11 

2.438E-11 

6.346E-12 

128 

4.996E-12 

7.749E-13 

7.352E-13 

2.583E-13 

256 

1.604E-13 

2.698E-14 

2.256E-14 

8.995E-15 

512 

5.080E-15 

8.862E-16 

6.985E-16 

2.954E-16 

CR 


4.93 

4.54 

5.05 

4.54 


interpolation methods here in our schemes. For high dimensional problems, 
this would results in the tensor Lagrange interpolation methods, which may 
be time consuming. Thus, we would suggest more feasible techniques such as 
the sparse grid interpolation, RBF interpolation, etc. This would be part of 
our future studies. 

• Rigorous stability and convergence analysis. This is also our ongoing project. 
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